#include "average.h"
#include "../sort/sort.h"

namespace roto
{ 
	namespace average
	{
		template< typename ET >
		ET rangelim(ET x, ET lim1, ET lim2){
			x = (x < lim1 ? lim1 : x);
			x = (x > lim2 ? lim2 : x);
			return x;
		}

		uchar* averaging(uchar* F, uint width, uint height, int m, int n, const int type, float Q)
		{
			float avr, mn, savr;
			int x1, y1, m1 = m, n1 = n;
			float* median;
			uchar* T = new uchar[width * height];
			
			if(m1 % 2 == 0)m1 ++;
			if(n1 % 2 == 0)m1 ++;
			m1 = (m1 - 1) / 2;
			n1 = (n1 - 1) / 2;
			mn = (float)(m * n);
			median = new float[(int)mn];

			switch(type){
				case 1:
					for (uint i = 0; i < height; ++ i)
						for (uint j = 0; j < width; ++ j)
						{
							long ij = i * width + j;
							avr = 0;
							for (int x = j - m1; x <= j + m1; ++ x)
								for (int y = i - n1; y <= i + n1; ++ y)
								{
									long xy = rangelim<float>((float)y, .0f, height - 1.0f) * width + rangelim<float>((float)x, .0f, width - 1.0f);
									avr += F[xy];
								}
							T[ij] = avr / mn;
						}
					break;
				case 2:
					for (uint i = 0; i < height; ++ i)
						for (uint j = 0; j < width; ++ j)
						{
							long ij = i * width + j;
							avr = 1;
							for (int x = j - m1; x <= j + m1; ++ x)
								for (int y = i - n1; y <= i + n1; ++ y)
								{
									long xy = rangelim<float>((float)y, .0f, height - 1.0f) * width + rangelim<float>((float)x, .0f, width - 1.0f);
									avr *= F[xy];
								}
							avr = exp(log(avr) / mn);
							T[ij] = avr;
						}
					break;
				case 3:
					for (uint i = 0; i < height; ++ i)
						for (uint j = 0; j < width; ++ j)
						{
							long ij = i * width + j;
							avr = 0;
							for (int x = j - m1; x <= j + m1; ++ x)
								for (int y = i - n1; y <= i + n1; ++ y)
								{
									long xy = rangelim<float>((float)y, .0f, height - 1.0f) * width + rangelim<float>((float)x, .0f, width - 1.0f);
									avr += 1.0f / F[xy];
								}
							avr = mn / avr; 
							T[ij] = avr;
						}
					break;
				case 4:
					for (uint i = 0; i < height; ++ i)
						for (uint j = 0; j < width; ++ j)
						{
							long ij = i * width + j;
							savr = avr = 0;
							for (int x = j - m1; x <= j + m1; ++ x)
								for (int y = i - n1; y <= i + n1; ++ y)
								{
									x1 = (int)rangelim<float>((float)x, .0f, width - 1.0f);
									y1 = (int)rangelim<float>((float)y, .0f, height - 1.0f);
									long xy1 = y1 * width + x1;
									avr += (F[xy1] != .0f ? ((Q + 1.0f) != 0 ? exp((Q + 1.0f) * log((double)F[xy1])) : 1.0f) : .0f);
									savr += (F[xy1] != .0f ? ((Q) != .0f ? exp((Q) * log((double)F[xy1])) : 1.0f) : 1.0f);
								}
							T[ij] = avr / savr;
						}
					break;
				case 5:
					for (uint i = 0; i < height; ++ i)
						for (uint j = 0; j < width; ++ j)
						{
							long ij = i * width + j;
							for (int x = j - m1; x <= j + m1; ++ x)
								for (int y = i - n1; y <= i + n1; ++ y)
								{
									long xy = rangelim<float>((float)y, .0f, height - 1.0f) * width + rangelim<float>((float)x, .0f, width - 1.0f);
									median[(x + m1 - i) * n  + y - j + n1] = F[xy];
								}
								sort::quickSortR<float>(median, mn - 1);
							if(Q == 1.0f)
								T[ij] = (median[0] + median[(uint)mn - 1]) * .5f;
							else
								T[ij] = median[(uint)((float)mn / 2.0f)];
							
						}
					break;
				default:
					return F;
			}
			return T;
		}
		
	}

}